%% tables_figures.m
%
% Create Tables and Figures for paper
%--------------------------------------------------------------------------

%%%%%%%%%%%%%%%%%%%%%%
saveFigs = 1;

benchmark_exp = 'beta1000rho125_outputAll';
benchmark_PB_P = 'beta830rho88_outputAll';

exp_proc = 'beta999rho133_outputAll';
exp_rational_slow = 'beta1000rho134_outputAll';
PB_noproc = 'beta675rho54_outputAll';
PB_rational_slow = 'beta725rho70_outputAll';
%%%%%%%%%%%%%%%%%%%%%%

%Set plot colors
colorB1 = [0 0 0];
colorPB = [0.6350 0.0780 0.1840];

addpath('subfunctions')


%--------------------------------------------------------------------------
%% Tables

    %Table of Steady State Moments (Tables 1 and 2)
    for runind = 1:2
        if runind == 1
            load(['output/', benchmark_exp]);
        elseif runind == 2
            load(['output/', benchmark_PB_P]);
        end
        g = g_stationary(:,:,:,rb_spot_ind_stationary);
        cStationary = cStoreSS{rb_spot_ind_stationary, rb_spot_ind_stationary};
        income_shares = sum(squeeze(sum(sum(g(:,:,:,:)*da*db))),2);

        avgLTV = sum(sum(sum(g.*(aaa./h)*da*db)));
        avgCC = sum(sum(sum(g(1:bzeroInd-1,:,:).*bbb(1:bzeroInd-1,:,:)*da*db)));
        shareCC = sum(sum(sum(g(1:bzeroInd-1,:,:)*da*db)));
        shareCCMax = sum(sum(g(1,:,:)*da*db));    

        durableShare = 0.125;
        depr = 0.22;
        rCurr = rbAll(rb_spot_ind_stationary);

        %Average quarterly MPC/MPX over small shock
        MPCshock = 1000/averageIncomeSCF; 
            consumption_FK;
            distShockType = 2;
            liqShock = MPCshock;
            distShockInterventionMatrix;
        MPC3_0pt25 = (SC*C_FK3_0pt25(:) - C_FK3_0pt25(:))./liqShock;
        MPC3_0pt25 = reshape(MPC3_0pt25, Nb,Na,Nz,Nr);
        avgMPC3_0pt25 = sum(sum(sum(sum(MPC3_0pt25.*g_stationary*da*db))));
            avgMPC3_0pt25_LMH = squeeze(sum(sum(sum(MPC3_0pt25.*g_stationary*da*db,4),2),1))./income_shares;
        Deriv_C3_0pt25 = (SC*Ec_FK3_0pt25(:) - Ec_FK3_0pt25(:))./liqShock;
        Deriv_C3_0pt25 = reshape(Deriv_C3_0pt25, Nb,Na,Nz,Nr);
        MPX3_0pt25 = (1-durableShare+((depr*durableShare)/(depr+rCurr)))*MPC3_0pt25 + (durableShare/(depr+rCurr))*Deriv_C3_0pt25;
        avgMPX3_0pt25 = sum(sum(sum(sum(MPX3_0pt25.*g_stationary*da*db))));

        %Average quarterly MPC/MPX over large shock        
        MPCshock = 10000/averageIncomeSCF; 
            distShockType = 2;
            liqShock = MPCshock;
            distShockInterventionMatrix;
        MPC3_0pt25 = (SC*C_FK3_0pt25(:) - C_FK3_0pt25(:))./liqShock;
        MPC3_0pt25 = reshape(MPC3_0pt25, Nb,Na,Nz,Nr);
        avgMPC3_0pt25_large = sum(sum(sum(sum(MPC3_0pt25.*g_stationary*da*db))));
        Deriv_C3_0pt25 = (SC*Ec_FK3_0pt25(:) - Ec_FK3_0pt25(:))./liqShock;
        Deriv_C3_0pt25 = reshape(Deriv_C3_0pt25, Nb,Na,Nz,Nr);
        MPX3_0pt25 = (1-durableShare+((depr*durableShare)/(depr+rCurr)))*MPC3_0pt25 + (durableShare/(depr+rCurr))*Deriv_C3_0pt25;
        avgMPX3_0pt25_large = sum(sum(sum(sum(MPX3_0pt25.*g_stationary*da*db))));

        preshock_table(:,runind) = [avgLTV; avgCC; avgMPC3_0pt25; avgMPC3_0pt25_LMH; avgMPC3_0pt25_large; avgMPX3_0pt25; avgMPX3_0pt25_large; shareCC; shareCCMax];
    end
    if saveFigs == 1
        save('figs/Tables/steady_state_table.mat', 'preshock_table')
    end

    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    
    %FP Consumption Response (Figure 4)
    FPshock_table = zeros(3,2);
    load(['output/', benchmark_exp]);
    for runind = 1:2
        if runind == 1
            load(['output/', benchmark_exp]);
            c_preshock = sum(meanc_t_store(:,1).*income_shares);
            load(['output/', benchmark_exp, '_FP_liquid']);
        elseif runind == 2
            load(['output/', benchmark_PB_P]);
            c_preshock = sum(meanc_t_store(:,1).*income_shares);
            load(['output/', benchmark_PB_P, '_FP_liquid']);
        end

        ind4Q = find(time_path >= 0.99,1)-1;
        ind8Q = find(time_path >= 1.99,1)-1;
        ind12Q = find(time_path >= 2.99,1)-1;

        Ccumulative = cumsum(Delta_path.*(sum(meanc_t_store(:,2:end).*income_shares)-c_preshock))/liqShock;
        FPshock_table(1,runind) = Ccumulative(ind4Q);
        FPshock_table(2,runind) = Ccumulative(ind8Q);
        FPshock_table(3,runind) = Ccumulative(ind12Q);
    end
    FPshock_table = round(FPshock_table,2);
    if saveFigs == 1
        save('figs/Tables/FP_table.mat', 'FPshock_table')
    end
    
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    %MP Present Value (Figure 5)
    MP_PV_table = zeros(3,2);
    for runind = 1:2
        if runind == 1
            load(['output/', benchmark_exp]);
        elseif runind == 2
            load(['output/', benchmark_PB_P]);
        end

        discountPath = exp(cumsum(-rbAll(rb_spot_ind_path).*Delta_path));
        integratedCPost = cumsum(sum(meanc_t_store(:,2:end).*income_shares).*Delta_path.*discountPath);
        integratedCPre = cumsum(sum(meanc_t_store(:,1).*income_shares).*Delta_path.*discountPath);
        CPreRate = income_shares' * meanc_t_store(:,1);

        ind4Q = find(time_path >= 0.99,1)-1;
        ind8Q = find(time_path >= 1.99,1)-1;
        ind12Q = find(time_path >= 2.99,1)-1;

        MP_PV_table(:,runind) = [(integratedCPost(ind4Q)-integratedCPre(ind4Q))./CPreRate; ...
                                (integratedCPost(ind8Q)-integratedCPre(ind8Q))./CPreRate; ...
                                (integratedCPost(ind12Q)-integratedCPre(ind12Q))./CPreRate];
        MP_PV_table(:,runind) = 100*MP_PV_table(:,runind);
    end
    if saveFigs == 1
        save('figs/Tables/MP_PV_table.mat', 'MP_PV_table')
    end

    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    
    %MP Decomposition (Table 3)
    cDecompositionTable = zeros(6,2);
    for runind = 1:2
        if runind == 1
            load(['output/', benchmark_exp]);
            cDecompositionTable(3,runind) = 100*(sum(meanc_t_store(:,2).*income_shares)./sum(meanc_t_store(:,1).*income_shares) - 1);
            load(['output/', benchmark_exp, '_MPdecomp1']);
            cDecompositionTable(1,runind) = 100*(sum(meanc_t_store(:,2).*income_shares)./sum(meanc_t_store(:,1).*income_shares) - 1);
            load(['output/', benchmark_exp, '_MPdecomp2']);
            cDecompositionTable(2,runind) = 100*(sum(meanc_t_store(:,2).*income_shares)./sum(meanc_t_store(:,1).*income_shares) - 1);
        elseif runind == 2
            load(['output/', benchmark_PB_P]);
            cDecompositionTable(3,runind) = 100*(sum(meanc_t_store(:,2).*income_shares)./sum(meanc_t_store(:,1).*income_shares) - 1);
            load(['output/', benchmark_PB_P, '_MPdecomp1']);
            cDecompositionTable(1,runind) = 100*(sum(meanc_t_store(:,2).*income_shares)./sum(meanc_t_store(:,1).*income_shares) - 1);
            load(['output/', benchmark_PB_P, '_MPdecomp2']);
            cDecompositionTable(2,runind) = 100*(sum(meanc_t_store(:,2).*income_shares)./sum(meanc_t_store(:,1).*income_shares) - 1);
        end
    end
    cDecompositionTable(4,:) = round(cDecompositionTable(1,:)./cDecompositionTable(3,:), 2);
    cDecompositionTable(5,:) = round(cDecompositionTable(2,:)./cDecompositionTable(3,:), 2);
    cDecompositionTable(6,:) = round(cDecompositionTable(3,:)./cDecompositionTable(3,:), 2);
    if saveFigs == 1
        save('figs/Tables/MP_decomposition_table.mat', 'cDecompositionTable')
    end
    
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    %MP Summary Moments (Appendix Table 6)
    for runind = 1:2
        if runind == 1
            load(['output/', benchmark_exp]);
        elseif runind == 2
            load(['output/', benchmark_PB_P]);
        end

        ind1Q = find(time_path >= 0.24,1)-1;
        ind4Q = find(time_path >= 0.99,1)-1;
        ind8Q = find(time_path >= 1.99,1)-1;
        ind12Q = find(time_path >= 2.99,1)-1;
        refi1Q = sum(sum(sum(g_t_store(:,:,:,rb_spot_ind_path(1),ind1Q)*da*db)));
        refi4Q = sum(sum(sum(g_t_store(:,:,:,rb_spot_ind_path(1),ind4Q)*da*db)));
        refi8Q = sum(sum(sum(g_t_store(:,:,:,rb_spot_ind_path(1),ind8Q)*da*db)));
        refi12Q = sum(sum(sum(g_t_store(:,:,:,rb_spot_ind_path(1),ind12Q)*da*db)));

        adj = adjStore{rb_spot_ind_stationary-1, rb_spot_ind_stationary};
        adjRegion = reshape(adj,Nb,Na,Nz).*ones(Nb,Na,Nz);
        aAdj = aAdj_combineStore{rb_spot_ind_stationary-1, rb_spot_ind_stationary};
            aAdj = reshape(aAdj,Nb,Na,Nz);
            %Cash-out defined as having bigger mortgage (with 5% cushion)
            cashoutFlag = reshape(aAdj,Nb,Na,Nz)>(1.05*aaa);
        bAdj = bAdj_combineStore{rb_spot_ind_stationary-1, rb_spot_ind_stationary};
            bAdj = reshape(bAdj,Nb,Na,Nz);

        refi = rewrite_combineStore{rb_spot_ind_stationary-1, rb_spot_ind_stationary};
            refi = adjRegion.*reshape(refi,Nb,Na,Nz).*ones(Nb,Na,Nz);

        shareRefiImpact = sum(sum(sum(refi.*g_stationary(:,:,:,rb_spot_ind_stationary)*da*db)));
        shareCashoutImpact = sum(sum(sum(cashoutFlag.*refi.*g_stationary(:,:,:,rb_spot_ind_stationary)*da*db)))./shareRefiImpact;

        bRefi = (bAdj-bbb).*refi;
        meanbRefi = (sum(sum(sum(bRefi.*g_stationary(:,:,:,rb_spot_ind_stationary)*da*db))))./shareRefiImpact;

        MPshock_table(:,runind) = [shareRefiImpact; shareCashoutImpact; refi1Q; refi4Q; refi8Q; refi12Q; meanbRefi];
    end
    if saveFigs == 1
        save('figs/Tables/MP_details_table.mat', 'MPshock_table')
    end
    
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    
%--------------------------------------------------------------------------
%% Figures -- Headline Figures

    %Monetary Policy (Figure 5)
    fig = figure; hold on;
        load(['output/', benchmark_exp]);
        plot(time_path(2:end), 100*(sum(meanc_t_store(:,2:end).*income_shares)./sum(meanc_t_store(:,1).*income_shares) - 1), 'LineWidth', 5, 'color', colorB1); hold on;
        
        load(['output/', benchmark_PB_P]);
        plot(time_path(2:end), 100*(sum(meanc_t_store(:,2:end).*income_shares)./sum(meanc_t_store(:,1).*income_shares) - 1), 'LineWidth', 5, 'color', colorPB); hold on;
        
        switch_one_q = 1;
        simulateForward;
        plot(time_path(2:end), 100*(sum(meanc_t_store(:,2:end).*income_shares)./sum(meanc_t_store(:,1).*income_shares) - 1), 'LineWidth', 5, 'color', [colorPB, 0.3]); hold on;
        
        legend('Exponential','Present Bias','Present Bias (1Q No Proc.)','FontSize', 20, 'Interpreter', 'Latex');
        xlim([0,3]); ylim([0, 7]); grid on;
        set(gca,'FontSize', 14);
        xlabel('Years','FontSize', 20, 'Interpreter', 'Latex'); ylabel('Consumption Elasticity','FontSize', 20, 'Interpreter', 'Latex');
        set(gcf,'PaperPosition',[0 0 8 6]); clear switch_one_q;
    if saveFigs == 1
        saveas(fig, ['figs/Policy_Functions/Figure5.png'])
    end
    
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    
    %Fiscal Policy (Figure 4)
    fig = figure; hold on; 
        load(['output/', benchmark_exp]);
        c_preshock = sum(meanc_t_store(:,1).*income_shares);
        load(['output/', benchmark_exp, '_FP_liquid']);
        plot(time_path(2:end), (sum(meanc_t_store(:,2:end).*income_shares)-c_preshock)/liqShock, 'LineWidth', 5, 'color', colorB1); hold on;
        
        load(['output/', benchmark_PB_P]);
        c_preshock = sum(meanc_t_store(:,1).*income_shares);
        load(['output/', benchmark_PB_P, '_FP_liquid']);
        plot(time_path(2:end), (sum(meanc_t_store(:,2:end).*income_shares)-c_preshock)/liqShock, 'LineWidth', 5, 'color', colorPB); hold on;
         
        xlim([0,3]); ylim([0,2]); grid on;
        set(gca,'FontSize', 14);
        xlabel('Years','FontSize', 20, 'Interpreter', 'Latex'); 
        ylabel('Consumption IRF','FontSize', 20, 'Interpreter', 'Latex');
        legend('Exponential','Present Bias','FontSize', 20, 'Interpreter', 'Latex');
        set(gcf,'PaperPosition',[0 0 8 6])
    if saveFigs == 1
        saveas(fig, ['figs/Policy_Functions/Figure4.png']);
    end
    
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    

%--------------------------------------------------------------------------
%% Figures -- Model Steady State Solution Figures

    %Consumption (Figure 2)     
    for figind = 1:2        
        if figind == 1
            load(['output/', benchmark_exp]);
            startPos = 1;
        elseif figind == 2
            load(['output/', benchmark_PB_P]);
            startPos = 2;
        end
        color0 = [0 0.24 0.64]; color40 = [0.8 0.3 0.08]; color80 = [0.6 0.81 0.36];
        cStat = cStoreSS{rb_spot_ind_stationary,rb_spot_ind_stationary};
            if slow_refi == 0
                adj = adjStore{rb_spot_ind_stationary, rb_spot_ind_stationary};
                adjRegion = reshape(adj,Nb,Na,Nz).*ones(Nb,Na,Nz);
                cStat(adjRegion == 1)=NaN;
            end
            if slow_refi == 1
                adj = adjStore{rb_spot_ind_stationary, rb_spot_ind_stationary};
                adjRegion = reshape(adj,Nb,Na,Nz).*ones(Nb,Na,Nz);
                cStatadj = cStat;
                cStat(adjRegion == 1)=NaN;
                cStatadj(adjRegion ~= 1)=NaN;
            end
        fig = figure;
        set(gcf,'PaperPosition',[0 0 12 5])
        subplot(1,3,1); hold on; grid on;
            plot(b(startPos:end),cStat(startPos:end,discretize(0, a./h),1), 'linewidth', 3, 'color', color0);
            plot(b(startPos:end),cStat(startPos:end,discretize(0.4, a./h),1), 'linewidth', 3, 'color', color40);
            plot(b(startPos:end),cStat(startPos:end,Na,1), 'linewidth', 3, 'color', color80); assert(a(Na)./h == 0.8);
            if slow_refi == 1
                plot(b(startPos:end),cStatadj(startPos:end,discretize(0, a./h),1), 'linewidth', 3, 'color', color0, 'linestyle', ':');
                plot(b(startPos:end),cStatadj(startPos:end,discretize(0.4, a./h),1), 'linewidth', 3, 'color', color40, 'linestyle', ':');
                plot(b(startPos:end),cStatadj(startPos:end,Na,1), 'linewidth', 3, 'color', color80, 'linestyle', ':');
            end
            plot(b(1),cStat(1,discretize(0, a./h),1), 'o', 'color', color0, 'markerfacecolor', color0, 'markersize', 12);
            plot(b(1),cStat(1,discretize(0.4, a./h),1), 'o', 'color', color40, 'markerfacecolor', color40, 'markersize', 12);
            plot(b(1),cStat(1,Na,1), 'o', 'color', color80, 'markerfacecolor', color80, 'markersize', 12);
            if slow_refi == 1
                plot(b(1),cStatadj(1,discretize(0, a./h),1), 'o', 'color', color0, 'markersize', 12);
                plot(b(1),cStatadj(1,discretize(0.4, a./h),1), 'o', 'color', color40, 'markersize', 12);
                plot(b(1),cStatadj(1,Na,1), 'o', 'color', color80, 'markersize', 12);
            end
            set(gca, 'FontSize', 12); ylim([.4, 1.3]);
            xlim([bmin,1.01]); legend('Mortgage (LTV) = 0','Mortgage (LTV) = 0.4','Mortgage (LTV) = 0.8', 'Interpreter', 'Latex', 'location', 'southeast', 'FontSize', 16);
            title('Low Income', 'Interpreter', 'Latex', 'FontSize', 18);
            xlabel('Liquid Wealth','FontSize',16, 'interpreter','latex');
            ylabel('Consumption','FontSize',16, 'interpreter','latex');
        subplot(1,3,2); hold on; grid on;
            plot(b(startPos:end),cStat(startPos:end,discretize(0, a./h),2), 'linewidth', 3, 'color', color0);
            plot(b(startPos:end),cStat(startPos:end,discretize(0.4, a./h),2), 'linewidth', 3, 'color', color40);
            plot(b(startPos:end),cStat(startPos:end,Na,2), 'linewidth', 3, 'color', color80);
            if slow_refi == 1
                plot(b(startPos:end),cStatadj(startPos:end,discretize(0, a./h),2), 'linewidth', 3, 'color', color0, 'linestyle', ':');
                plot(b(startPos:end),cStatadj(startPos:end,discretize(0.4, a./h),2), 'linewidth', 3, 'color', color40, 'linestyle', ':');
                plot(b(startPos:end),cStatadj(startPos:end,Na,2), 'linewidth', 3, 'color', color80, 'linestyle', ':');
            end
            plot(b(1),cStat(1,discretize(0, a./h),2), 'o', 'color', color0, 'markerfacecolor', color0, 'markersize', 12);
            plot(b(1),cStat(1,discretize(0.4, a./h),2), 'o', 'color', color40, 'markerfacecolor', color40, 'markersize', 12);
            plot(b(1),cStat(1,Na,2), 'o', 'color', color80, 'markerfacecolor', color80, 'markersize', 12);
            if slow_refi == 1
                plot(b(1),cStatadj(1,discretize(0, a./h),2), 'o', 'color', color0, 'markersize', 12);
                plot(b(1),cStatadj(1,discretize(0.4, a./h),2), 'o', 'color', color40, 'markersize', 12);
                plot(b(1),cStatadj(1,Na,2), 'o', 'color', color80, 'markersize', 12);
            end
            set(gca, 'FontSize', 12); ylim([.4, 1.3]);
            xlim([bmin,1.01]); legend('Mortgage (LTV) = 0','Mortgage (LTV) = 0.4','Mortgage (LTV) = 0.8', 'Interpreter', 'Latex', 'location', 'southeast', 'FontSize', 16); 
            title('Middle Income', 'Interpreter', 'Latex', 'FontSize', 18);
            xlabel('Liquid Wealth','FontSize',16, 'interpreter','latex');
            ylabel('Consumption','FontSize',16, 'interpreter','latex');
        subplot(1,3,3); hold on; grid on;
            plot(b(startPos:end),cStat(startPos:end,discretize(0, a./h),3), 'linewidth', 3, 'color', color0);
            plot(b(startPos:end),cStat(startPos:end,discretize(0.4, a./h),3), 'linewidth', 3, 'color', color40);
            plot(b(startPos:end),cStat(startPos:end,Na,3), 'linewidth', 3, 'color', color80);
            if slow_refi == 1
                plot(b(startPos:end),cStatadj(startPos:end,discretize(0, a./h),3), 'linewidth', 3, 'color', color0, 'linestyle', ':');
                plot(b(startPos:end),cStatadj(startPos:end,discretize(0.4, a./h),3), 'linewidth', 3, 'color', color40, 'linestyle', ':');
                plot(b(startPos:end),cStatadj(startPos:end,Na,3), 'linewidth', 3, 'color', color80, 'linestyle', ':');
            end
            plot(b(1),cStat(1,discretize(0, a./h),3), 'o', 'color', color0, 'markerfacecolor', color0, 'markersize', 12);
            plot(b(1),cStat(1,discretize(0.4, a./h),3), 'o', 'color', color40, 'markerfacecolor', color40, 'markersize', 12);
            plot(b(1),cStat(1,Na,3), 'o', 'color', color80, 'markerfacecolor', color80, 'markersize', 12);
            if slow_refi == 1
                plot(b(1),cStatadj(1,discretize(0, a./h),3), 'o', 'color', color0, 'markersize', 12);
                plot(b(1),cStatadj(1,discretize(0.4, a./h),3), 'o', 'color', color40, 'markersize', 12);
                plot(b(1),cStatadj(1,Na,3), 'o', 'color', color80, 'markersize', 12);
            end
            set(gca, 'FontSize', 12); ylim([.4, 1.3]);
            xlim([bmin,1.01]); legend('Mortgage (LTV) = 0','Mortgage (LTV) = 0.4','Mortgage (LTV) = 0.8', 'Interpreter', 'Latex', 'location', 'southeast', 'FontSize', 16);
            title('High Income', 'Interpreter', 'Latex', 'FontSize', 18);
            xlabel('Liquid Wealth','FontSize',16, 'interpreter','latex');
            ylabel('Consumption','FontSize',16, 'interpreter','latex');
        if saveFigs == 1
            set(gcf,'PaperPosition',[0 0 16 5])
            saveas(fig, ['figs/Policy_Functions/beta', num2str(1000*beta),  '_Figure2.png'])
        end
    end
    
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    
    %Overall Distribution (Appendix Figure 14)
    for figind = 1:2
        if figind == 1
            load(['output/', benchmark_exp]);
        elseif figind == 2
            load(['output/', benchmark_PB_P]);
        end
        fig=figure; set(gcf,'PaperPosition',[0 0 16 5])
        g = g_stationary;
        for z_ind = 1:Nz
            subplot(1,Nz,z_ind)
            surf(b,a./h,g(:, :, z_ind,rb_spot_ind_stationary)'./income_shares(z_ind), 'LineStyle', ':')
            view([0 90])
            caxis([0 1.25]); 
            xlabel('Liquid Wealth','FontSize',16, 'interpreter','latex')
            ylabel('Mortgage (LTV)','FontSize',16, 'interpreter','latex')
            xlim([bmin, 1]); ylim([amin, amax./h]); 
            if z_ind == 1
                title('Low Income','FontSize',17, 'interpreter','latex');
            elseif z_ind == 2
                title('Middle Income','FontSize',17, 'interpreter','latex');
            elseif z_ind == 3
                title('High Income','FontSize',17, 'interpreter','latex');
            end
        end
        if saveFigs == 1
            saveas(fig, ['figs/Distributions/beta', num2str(1000*beta),  '_AFigure14.png'])
        end
    end
    
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    
    %Liquid Wealth Distribution (Figure 3)
    fig=figure; hold on;
    for figind = 1:2
        if figind == 1
            load(['output/', benchmark_exp]);
        elseif figind == 2
            load(['output/', benchmark_PB_P]);
        end
        
        dist_b = sum(sum(g_stationary(:,:,:,rb_spot_ind_stationary), 3), 2);
        
        if figind == 1
            plot(b, dist_b, 'Linewidth', 5, 'color', colorB1);
        elseif figind == 2
            plot(b,dist_b, 'Linewidth', 5, 'color', colorPB);
            shareMinP = round(dist_b(1)*da*db*100, 1);
        end
    end
    ymax = shareMinP/(da*db*300);
    xlim([bmin,1]); ylim([0, ymax]); grid on; 
    set(gca,'FontSize',12); set(gca, 'YTickLabel', []);
    xlabel('Liquid Wealth','FontSize', 20, 'Interpreter', 'Latex'); ylabel('Marginal Distribution','FontSize', 20, 'Interpreter', 'Latex');
    legend('Exponential','Present Bias','FontSize', 20, 'Interpreter', 'Latex');
    if saveFigs == 1
        set(gcf,'PaperPosition',[0 0 8 6])
        saveas(fig, ['figs/Distributions/Figure3.png']);
    end
   
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    
    %Mortgage Distribution (Appendix Figure 13)
    fig=figure; hold on;
    set(gcf,'PaperPosition',[0 0 9 6])
    for figind = 1:2
        if figind == 1
            load(['output/', benchmark_exp]);
        elseif figind == 2
            load(['output/', benchmark_PB_P]);
        end
        
        dist_a = squeeze(sum(sum(g_stationary(:,:,:,rb_spot_ind_stationary), 3), 1));
        
        if figind == 1
            plot(a./h, dist_a, 'Linewidth', 5, 'color', colorB1); hold on; grid on;
        elseif figind == 2
            plot(a./h,dist_a, 'Linewidth', 5, 'color', colorPB);
        end
    end
    xlim([0,thetam]); ylim([0, 100]); grid on; 
    set(gca,'FontSize',12); set(gca, 'YTickLabel', []);
    xlabel('Mortgage (LTV)','FontSize', 20, 'Interpreter', 'Latex'); ylabel('Marginal Distribution','FontSize', 20, 'Interpreter', 'Latex');
    legend('Exponential','Present Bias','FontSize', 20, 'Interpreter', 'Latex');
    if saveFigs == 1
        saveas(fig, ['figs/Distributions/AFigure13.png'])
    end
    
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    
    %MPCs and Liquidity (Figure 3)
    MPC3_0pt25_b = zeros(Nb,2);
    for figind = 1:2
        if figind == 1
            load(['output/', benchmark_exp]);
        elseif figind == 2
            load(['output/', benchmark_PB_P]);
        end
        MPCshock = 1000/averageIncomeSCF; 
        consumption_FK;
        distShockType = 2;
        liqShock = MPCshock;
        distShockInterventionMatrix;
        MPC3_0pt25 = (SC*C_FK3_0pt25(:) - C_FK3_0pt25(:))./liqShock;
        MPC3_0pt25 = reshape(MPC3_0pt25, Nb,Na,Nz,Nr);
        MPC3_0pt25_b(:,figind) = (sum(sum((MPC3_0pt25(:,:,:,rb_spot_ind_stationary).*g_stationary(:,:,:,rb_spot_ind_stationary)),3),2))./sum(sum(g_stationary(:,:,:,rb_spot_ind_stationary),3),2);
    end
    fig = figure; hold on; grid on;
        plot(b, MPC3_0pt25_b(:,1), 'Linewidth', 5, 'color', colorB1); 
        plot(b, MPC3_0pt25_b(:,2), 'Linewidth', 5, 'color', colorPB);
        set(gca,'FontSize',12);
        xlim([bmin,1]); ylim([0, 0.7]);
        xlabel('Liquid Wealth','FontSize', 20, 'Interpreter', 'Latex'); ylabel('MPC','FontSize', 20, 'Interpreter', 'Latex');
        legend('Exponential','Present Bias','FontSize', 20, 'Interpreter', 'Latex');
        set(gcf,'PaperPosition',[0 0 8 6])
    if saveFigs == 1
        saveas(fig, ['figs/Policy_Functions/Figure3.png']);
    end
    
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    %MPCs by shock size (Appendix Figure 12)
    fig=figure; hold on;
    shockLevels = [1000, 5000:5000:50000]./averageIncomeSCF;
    MPC_size_table = zeros(length(shockLevels), 2);
    for runind = 1:2
        if runind == 1
            load(['output/', benchmark_exp]);
        elseif runind == 2
            load(['output/', benchmark_PB_P]);
        end
        consumption_FK;
    
        for sL = 1:length(shockLevels)
            distShockType = 2;
            liqShock = shockLevels(sL);
            distShockInterventionMatrix;    
            MPC3_0pt25 = (SC*C_FK3_0pt25(:) - C_FK3_0pt25(:))./liqShock;
            MPC_size_table(sL,runind) = sum(sum(sum(sum(g_stationary.*reshape(MPC3_0pt25, Nb,Na,Nz,Nr)*da*db))));
        end
        if runind == 1
            plot(shockLevels*averageIncomeSCF, MPC_size_table(:,runind),'-o', 'Linewidth', 3, 'color', colorB1, 'MarkerSize',11); hold on; grid on;
        elseif runind == 2
            plot(shockLevels*averageIncomeSCF, MPC_size_table(:,runind), '-o','Linewidth', 3, 'color', colorPB, 'MarkerSize',11);
        end 
    end
    set(gca,'FontSize',14);
    xlim([0,50000]); ylim([0, 0.15]);
    xlabel('Transfer Amount','FontSize', 20, 'Interpreter', 'Latex'); ylabel('Quarterly MPC','FontSize', 20, 'Interpreter', 'Latex');
    legend('Exponential','Present Bias','FontSize', 20, 'Interpreter', 'Latex');
    xtickformat('$%,.0f');
    ax = gca;
    ax.XAxis.Exponent = 0;
    set(gcf,'PaperPosition',[0 0 9 6])
    if saveFigs == 1
        saveas(fig, ['figs/Policy_Functions/AFigure12_1.png'])
    end
    
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    
    %MPCs by shock size (over b) (Appendix Figure 12)
    fig=figure; hold on;
    shockLevels = [1000/averageIncomeSCF, 10000/averageIncomeSCF, 25000/averageIncomeSCF];
    load(['output/', benchmark_PB_P]);
    consumption_FK;
    for sL = 1:length(shockLevels) 
        MPCshock = shockLevels(sL); 
        distShockType = 2;
        liqShock = MPCshock;
        distShockInterventionMatrix;
        MPC3_0pt25 = (SC*C_FK3_0pt25(:) - C_FK3_0pt25(:))./liqShock;
        MPC3_0pt25 = reshape(MPC3_0pt25, Nb,Na,Nz,Nr);
        MPC3_0pt25_b = (sum(sum((MPC3_0pt25(:,:,:,rb_spot_ind_stationary).*g_stationary(:,:,:,rb_spot_ind_stationary)),3),2))./sum(sum(g_stationary(:,:,:,rb_spot_ind_stationary),3),2);
        if sL == 1
            plot(b, MPC3_0pt25_b, 'Linewidth', 5, 'color', colorPB); hold on; grid on;
        elseif sL == 2
            plot(b, MPC3_0pt25_b, 'Linewidth', 3, 'color', [.5, .5, .5]);
        else
            plot(b, MPC3_0pt25_b, 'Linewidth', 3, 'color', [.5, .5, .5], 'linestyle', ':');
        end
    end
    set(gca,'FontSize',12);
    xlim([bmin,0.5]); ylim([0, 0.7]);
    xlabel('Liquid Wealth','FontSize', 20, 'Interpreter', 'Latex'); ylabel('MPC','FontSize', 20, 'Interpreter', 'Latex');
    legend('Transfer $= \$1,000$','Transfer $= \$10,000$','Transfer $= \$25,000$','FontSize', 20, 'Interpreter', 'Latex');
    set(gcf,'PaperPosition',[0 0 8 6])
    if saveFigs == 1
        saveas(fig, ['figs/Policy_Functions/AFigure12_2.png']);
    end
    
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    
    %Phase Diagrams (Figure 1 and Appendix Figure 16)
    for figind = 1:2
        if figind == 1
            load(['output/', benchmark_exp]);
        elseif figind == 2
            load(['output/', benchmark_PB_P]);
        end
        
        for rb_ind = rb_spot_ind_stationary-1:rb_spot_ind_stationary
            fig=figure;
                adj = adjStore{rb_ind, rb_spot_ind_stationary};
                aAdj = aAdj_combineStore{rb_ind, rb_spot_ind_stationary};
                bAdj = bAdj_combineStore{rb_ind, rb_spot_ind_stationary};
                rewrite_combine = rewrite_combineStore{rb_ind,rb_spot_ind_stationary};
                cashoutFlag = reshape(aAdj,Nb,Na,Nz)>(1.05*aaa);
                rewriteFlag = reshape(rewrite_combine,Nb,Na,Nz);
                rateFlag = logical(rewriteFlag.*(1-cashoutFlag));
                adjRegion = reshape(adj,Nb,Na,Nz).*ones(Nb,Na,Nz);
                adjRegionCO = logical(adjRegion.*cashoutFlag);
                adjRegionRate = logical(adjRegion.*rateFlag);
                adjRegionPrepay = logical(adjRegion.*(1-rewriteFlag));        

                sStat = sssStoreSS{rb_ind,rb_spot_ind_stationary};
                if slow_refi == 0
                    sStat(adjRegion == 1)=NaN;
                end
                if slow_refi == 1
                    sStatadj = sStat;
                    sStat(adjRegion == 1)=NaN;
                    sStatadj(adjRegion ~= 1)=NaN;
                end
                
                bb = bbb(:,:,1);
                aa = aaa(:,:,1);    
                bb = bb([1:20:Nb],[1:5:Na,Na]);
                aa = aa([1:20:Nb],[1:5:Na,Na]);
                sStat = sStat([1:20:Nb],[1:5:Na,Na],:);
                if slow_refi == 1
                    sStatadj = sStatadj([1:20:Nb],[1:5:Na,Na],:);
                end

                for z_ind = 1:Nz
                    subplot(1,Nz,z_ind)
                        %cashout = red; rate refi = gray; prepay = blue
                        CO = ones(Nb,Na,3);
                        CO(:,:,1) = ones(Nb,Na)*0.95;
                        CO(:,:,2) = ones(Nb,Na)*0.95;
                        CO(:,:,3) = ones(Nb,Na)*0.95;
                        CO(:,:,1) = CO(:,:,1) + adjRegionCO(:,:,z_ind)*(.5-.95);
                        CO(:,:,2) = CO(:,:,2) + adjRegionCO(:,:,z_ind)*(.03-.95);
                        CO(:,:,3) = CO(:,:,3) + adjRegionCO(:,:,z_ind)*(.1-.95);
                        CO(:,:,1) = CO(:,:,1) + adjRegionRate(:,:,z_ind)*(.4-.95);
                        CO(:,:,2) = CO(:,:,2) + adjRegionRate(:,:,z_ind)*(.4-.95);
                        CO(:,:,3) = CO(:,:,3) + adjRegionRate(:,:,z_ind)*(.4-.95);
                        CO(:,:,1) = CO(:,:,1) + adjRegionPrepay(:,:,z_ind)*(.62-.95);
                        CO(:,:,2) = CO(:,:,2) + adjRegionPrepay(:,:,z_ind)*(.84-.95);
                        CO(:,:,3) = CO(:,:,3) + adjRegionPrepay(:,:,z_ind)*(.93-.95);

                        surf(b,a./h,-1*adjRegion(:, :, z_ind)', permute(CO,[2,1,3]), 'FaceAlpha', 0.5); hold on;
                        view([0 90])

                        quiver(bb(:),aa(:)./h,reshape(sStat(:,:,z_ind), numel(aa),1),-xi*(aa(:)./h), 'linewidth', 1.5, 'Marker', '.', 'MarkerSize', 15);
                        if figind == 2
                           hold on;
                           quiver(bb(:),aa(:)./h,reshape(sStatadj(:,:,z_ind), numel(aa),1),-xi*(aa(:)./h), 'linewidth', 1, 'color', [.5 .5 .5], 'Marker', '.', 'MarkerSize', 15);
                        end

                        xlabel('Liquid Wealth','FontSize',16, 'interpreter','latex')
                        ylabel('Mortgage (LTV)','FontSize',16, 'interpreter','latex')
                        if z_ind == 1
                            title('Low Income','FontSize',18, 'interpreter','latex');
                        elseif z_ind == 2
                            title('Middle Income','FontSize',18, 'interpreter','latex');
                        elseif z_ind == 3
                            title('High Income','FontSize',18, 'interpreter','latex');
                        end
                        xlim([bmin, 1]); ylim([amin, amax./h]);
                        shading flat
                end
            if saveFigs == 1
                set(gcf,'PaperPosition',[0 0 14 4])
                saveas(fig, ['figs/Policy_Functions/beta', num2str(1000*beta), '_r', num2str(rb_ind),  '_Figure1_AFigure16.png'])
            end
        end
    end
    
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    
    
%--------------------------------------------------------------------------
%% Figures -- Extensions  
   
    %Distributional Effects (Appendix Figure 6)
    for runind = 1:2
        fig = figure; hold on; 
        if runind == 1
            load(['output/', benchmark_exp]);
            cc = colorB1;
        else
            load(['output/', benchmark_PB_P]);
            cc = colorPB;
        end
        c_pre_1 = cStore{rb_spot_ind_stationary,1};
        c_pre_2 = cStore{rb_spot_ind_stationary,2};
        c_pre_3 = cStore{rb_spot_ind_stationary,3};
        c_pre_4 = cStore{rb_spot_ind_stationary,4};
        c_pre_all = [c_pre_1(:); c_pre_2(:); c_pre_3(:); c_pre_4(:)];
        c_pre_all = reshape(c_pre_all,Nb,Na,Nz,Nr);
        
        distShockType = 2;
        distShockInterventionMatrix;
        c_post_all = SC*c_pre_all(:);
        c_post_all = reshape(c_post_all,Nb,Na,Nz,Nr);
        
        c_shock = (c_post_all(:,:,:,rb_spot_ind_stationary) - c_pre_all(:,:,:,rb_spot_ind_stationary))./liqShock;
        
        hstep = 0.025;
        cStat = cStoreSS{rb_spot_ind_stationary,rb_spot_ind_stationary};
        c_bins = [min(cStat(:))-hstep:hstep:max(cStat(:))+hstep];
        
        c_mass = zeros(length(c_bins)-1,1);
        deltaC = zeros(length(c_bins)-1,1);

        cboomind = discretize(cStat, c_bins);
        grel = g_stationary(:,:,:,rb_spot_ind_stationary);
        for ind = 1:Nb*Na*Nz
            indtmp = cboomind(ind);
            c_mass(indtmp) = c_mass(indtmp) + grel(ind);
            deltaC(indtmp) = deltaC(indtmp) + c_shock(ind)*grel(ind);
        end
        plot(c_bins(1:end-1), deltaC./c_mass, 'linewidth', 5.5, 'color', cc); hold on; grid on;
        set(gca,'FontSize', 14); 
        xlabel('Pre-Shock Consumption','FontSize', 20, 'Interpreter', 'Latex'); ylim([0,30]);
        ylabel('Consumption IRF','FontSize', 20, 'Interpreter', 'Latex');
        set(gcf,'PaperPosition',[0 0 8 6]);
        xlim([c_bins(find(c_mass > 0,1))-hstep/2,1.2]);
        yyaxis right;
        barplot = bar(c_bins(1:end-1), c_mass*da*db, 'FaceColor', cc);
            barplot.FaceAlpha = 0.1;
        ax = gca;
        ax.YColor = 'black';
        ylim([0,0.2]);
        legend('Consumption IRF', 'Distribution','FontSize', 18, 'Interpreter', 'Latex');
        if saveFigs == 1
            saveas(fig, ['figs/Policy_Functions/beta', num2str(1000*beta),  '_AFigure6_1.png'])
        end
    end
    
    for runind = 1:2
        fig = figure; hold on; 
        if runind == 1
            load(['output/', benchmark_exp]);
            cc = colorB1;
        else
            load(['output/', benchmark_PB_P]);
            cc = colorPB;
        end
        c_pre = cStore{rb_spot_ind_stationary, rb_spot_ind_stationary};
        c_post_1 = cStore{rb_spot_ind_stationary-1,1};
        c_post_2 = cStore{rb_spot_ind_stationary-1,2};
        c_post_3 = cStore{rb_spot_ind_stationary-1,3};
        c_post_4 = cStore{rb_spot_ind_stationary-1,4};

        %Get post-shock consumption for those who adjust due to the shock
        %(use intervention matrix C for rb_spot_ind = 2)
        c_post_all = [c_post_1(:); c_post_2(:); c_post_3(:); c_post_4(:)];
        c_plot_post = interventionStore{rb_spot_ind_stationary-1}*c_post_all;
            c_plot_post = reshape(c_plot_post,Nb,Na,Nz,Nr);
            %c_plot_post gives consumption just after shock
            c_plot_post = c_plot_post(:,:,:,rb_spot_ind_stationary);
        c_shock = 100*((c_plot_post./c_pre)-1);
        if slow_refi == 1
            %No adjustment on impact
            c_plot_post = reshape(c_post_all,Nb,Na,Nz,Nr);
            c_plot_post = c_plot_post(:,:,:,rb_spot_ind_stationary);
            c_shock = 100*((c_plot_post./c_pre)-1);
        end
        
        hstep = 0.025;
        cStat = cStoreSS{rb_spot_ind_stationary,rb_spot_ind_stationary};
        c_bins = [min(cStat(:))-hstep:hstep:max(cStat(:))+hstep];
        
        c_mass = zeros(length(c_bins)-1,1);
        deltaC = zeros(length(c_bins)-1,1);

        cboomind = discretize(cStat, c_bins);
        grel = g_stationary(:,:,:,rb_spot_ind_stationary);
        for ind = 1:Nb*Na*Nz
            indtmp = cboomind(ind);
            c_mass(indtmp) = c_mass(indtmp) + grel(ind);
            deltaC(indtmp) = deltaC(indtmp) + c_shock(ind)*grel(ind);
        end
        plot(c_bins(1:end-1), deltaC./c_mass, 'linewidth', 5.5, 'color', cc); hold on; grid on;
        set(gca,'FontSize', 14); 
        xlabel('Pre-Shock Consumption','FontSize', 20, 'Interpreter', 'Latex'); ylim([0,20]);
        ylabel('Consumption Elasticity','FontSize', 20, 'Interpreter', 'Latex');
        set(gcf,'PaperPosition',[0 0 8 6]);
        xlim([c_bins(find(c_mass > 0,1))-hstep/2,1.2]);
        yyaxis right;
        barplot = bar(c_bins(1:end-1), c_mass*da*db, 'FaceColor', cc);
            barplot.FaceAlpha = 0.1;
        ax = gca;
        ax.YColor = 'black';
        ylim([0,0.2]);
        legend('Consumption Elasticity', 'Distribution','FontSize', 18, 'Interpreter', 'Latex');
        if saveFigs == 1
            saveas(fig, ['figs/Policy_Functions/beta', num2str(1000*beta),  '_AFigure6_2.png'])
        end
    end
    
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
        
    %FP Comparison Liquid vs. Illiquid (Appendix Figure 15)
    fig = figure; hold on;
        load(['output/', benchmark_exp]);
        c_preshock = sum(meanc_t_store(:,1).*income_shares);
        load(['output/', benchmark_exp, '_FP_liquid']);
        plot(time_path(2:end), (sum(meanc_t_store(:,2:end).*income_shares)-c_preshock)/liqShock, 'LineWidth', 5, 'color', colorB1); hold on;
        
        load(['output/', benchmark_PB_P]);
        c_preshock = sum(meanc_t_store(:,1).*income_shares);
        load(['output/', benchmark_PB_P, '_FP_liquid']);
        plot(time_path(2:end), (sum(meanc_t_store(:,2:end).*income_shares)-c_preshock)/liqShock, 'LineWidth', 5, 'color', colorPB); hold on;
        
        xlim([0,3]); ylim([0,2]); grid on;
        set(gca,'FontSize', 14);
        xlabel('Years','FontSize', 20, 'Interpreter', 'Latex'); 
        ylabel('Consumption IRF','FontSize', 20, 'Interpreter', 'Latex');
        legend('Exponential','Present Bias','FontSize', 20, 'Interpreter', 'Latex');
        set(gcf,'PaperPosition',[0 0 6 6]);
    if saveFigs == 1
        saveas(fig, ['figs/Policy_Functions/AFigure15_1.png'])
    end
    
    fig = figure; hold on;
        load(['output/', benchmark_exp]);
        c_preshock = sum(meanc_t_store(:,1).*income_shares);
            notFullBenefit = find(a < -mortgageShock);
            shockAdjust = 0;
            for ind = 1:length(notFullBenefit)
                shockAdjust = shockAdjust + (-mortgageShock - a(ind))*(sum(sum(sum(g_stationary(:,ind,:,:)*da*db))));
            end
        load(['output/', benchmark_exp, '_FP_illiquid']);
        plot(time_path(2:end), (sum(meanc_t_store(:,2:end).*income_shares)-c_preshock)/(-mortgageShock-shockAdjust), 'LineWidth', 5, 'color', colorB1); hold on;
        
        load(['output/', benchmark_PB_P]);
        c_preshock = sum(meanc_t_store(:,1).*income_shares);
            notFullBenefit = find(a < -mortgageShock);
            shockAdjust = 0;
            for ind = 1:length(notFullBenefit)
                shockAdjust = shockAdjust + (-mortgageShock - a(ind))*(sum(sum(sum(g_stationary(:,ind,:,:)*da*db))));
            end
        load(['output/', benchmark_PB_P, '_FP_illiquid']);
        plot(time_path(2:end), (sum(meanc_t_store(:,2:end).*income_shares)-c_preshock)/(-mortgageShock-shockAdjust), 'LineWidth', 5, 'color', colorPB); hold on;
        
        xlim([0,3]); ylim([0,2]); grid on;
        set(gca,'FontSize', 14);
        xlabel('Years','FontSize', 20, 'Interpreter', 'Latex'); 
        ylabel('Consumption IRF','FontSize', 20, 'Interpreter', 'Latex');
        legend('Exponential','Present Bias','FontSize', 20, 'Interpreter', 'Latex');
        set(gcf,'PaperPosition',[0 0 6 6]);
    if saveFigs == 1
        saveas(fig, ['figs/Policy_Functions/AFigure15_2.png'])
    end
  
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
            
    %ARMs (Appendix Figure 17)
    fig = figure; hold on;
        load(['output/', benchmark_PB_P]);
        plot(time_path(2:end), 100*(sum(meanc_t_store(:,2:end).*income_shares)./sum(meanc_t_store(:,1).*income_shares) - 1), 'LineWidth', 5, 'color', colorPB); hold on;
        load(['output/', benchmark_PB_P, '_ARM']);
        plot(time_path(2:end), 100*(sum(meanc_t_store(:,2:end).*income_shares)./sum(meanc_t_store(:,1).*income_shares) - 1), 'LineWidth', 3, 'color', colorPB, 'linestyle', ':'); hold on;
        legend('Present Bias', 'Present Bias (ARM)', 'FontSize', 20, 'Interpreter', 'Latex');
        xlim([0,3]); ylim([0, 7]); grid on;
        set(gca,'FontSize', 14);
        xlabel('Years','FontSize', 20, 'Interpreter', 'Latex'); ylabel('Percent Change in Consumption','FontSize', 20, 'Interpreter', 'Latex');
        set(gcf,'PaperPosition',[0 0 8 6])
    if saveFigs == 1
        saveas(fig, ['figs/Policy_Functions/AFigure17.png']);
    end
    
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    
    %Monetary Policy with House Price Shocks (Appendix Figure 7)
    fig = figure; hold on;
        load(['output/', benchmark_exp, '_MP_HP_shock_25']);
        meanc_t_store_shock = meanc_t_store;
        load(['output/', benchmark_exp, '_HP_shock_25']);
        plot(time_path(2:end), 100*(sum(meanc_t_store_shock(:,2:end).*income_shares)./sum(meanc_t_store(:,2:end).*income_shares) - 1), 'LineWidth', 5, 'color', colorB1); hold on;
        load(['output/', benchmark_PB_P, '_MP_HP_shock_25']);
        meanc_t_store_shock = meanc_t_store;
        load(['output/', benchmark_PB_P, '_HP_shock_25']);
        plot(time_path(2:end), 100*(sum(meanc_t_store_shock(:,2:end).*income_shares)./sum(meanc_t_store(:,2:end).*income_shares) - 1), 'LineWidth', 5, 'color', colorPB); hold on;
        
        load(['output/', benchmark_exp]);
        tmp1=plot(time_path(2:end), 100*(sum(meanc_t_store(:,2:end).*income_shares)./sum(meanc_t_store(:,1).*income_shares) - 1), 'LineWidth', 3, 'color', colorB1); hold on;
        tmp1.Color(4) = 0.3;
        load(['output/', benchmark_PB_P]);
        tmp1=plot(time_path(2:end), 100*(sum(meanc_t_store(:,2:end).*income_shares)./sum(meanc_t_store(:,1).*income_shares) - 1), 'LineWidth', 3, 'color', colorPB); hold on;
        tmp1.Color(4) = 0.3;
        legend('Exponential','Present Bias','FontSize', 20, 'Interpreter', 'Latex');
        xlim([0,3]); ylim([0, 7]); grid on;
        set(gca,'FontSize', 14);
        xlabel('Years','FontSize', 20, 'Interpreter', 'Latex'); ylabel('Consumption Elasticity','FontSize', 20, 'Interpreter', 'Latex');
        set(gcf,'PaperPosition',[0 0 8 6])
    if saveFigs == 1
        saveas(fig, ['figs/Policy_Functions/AFigure7_1.png'])
    end
    
    fig = figure; hold on;
        load(['output/', benchmark_exp, '_MP_HP_shock_plus25']);
        meanc_t_store_shock = meanc_t_store;
        load(['output/', benchmark_exp, '_HP_shock_plus25']);
        plot(time_path(2:end), 100*(sum(meanc_t_store_shock(:,2:end).*income_shares)./sum(meanc_t_store(:,2:end).*income_shares) - 1), 'LineWidth', 5, 'color', colorB1); hold on;
        load(['output/', benchmark_PB_P, '_MP_HP_shock_plus25']);
        meanc_t_store_shock = meanc_t_store;
        load(['output/', benchmark_PB_P, '_HP_shock_plus25']);
        plot(time_path(2:end), 100*(sum(meanc_t_store_shock(:,2:end).*income_shares)./sum(meanc_t_store(:,2:end).*income_shares) - 1), 'LineWidth', 5, 'color', colorPB); hold on;
      
        load(['output/', benchmark_exp]);
        tmp1=plot(time_path(2:end), 100*(sum(meanc_t_store(:,2:end).*income_shares)./sum(meanc_t_store(:,1).*income_shares) - 1), 'LineWidth', 3, 'color', colorB1); hold on;
        tmp1.Color(4) = 0.3;
        load(['output/', benchmark_PB_P]);
        tmp1=plot(time_path(2:end), 100*(sum(meanc_t_store(:,2:end).*income_shares)./sum(meanc_t_store(:,1).*income_shares) - 1), 'LineWidth', 3, 'color', colorPB); hold on;
        tmp1.Color(4) = 0.3;
        legend('Exponential','Present Bias','FontSize', 20, 'Interpreter', 'Latex');
        xlim([0,3]); ylim([0, 7]); grid on;
        set(gca,'FontSize', 14);
        xlabel('Years','FontSize', 20, 'Interpreter', 'Latex'); ylabel('Consumption Elasticity','FontSize', 20, 'Interpreter', 'Latex');
        set(gcf,'PaperPosition',[0 0 8 6])
    if saveFigs == 1
        saveas(fig, ['figs/Policy_Functions/AFigure7_2.png'])
    end
    
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    
    %Fiscal Policy with House Price Shocks (Appendix Figure 8)
    fig = figure; hold on;     
        load(['output/', benchmark_exp, '_HP_shock_25']);
        meanc_t_store_pre = meanc_t_store;
        load(['output/', benchmark_exp, '_FP_liquid_HP_shock_25']);
        plot(time_path(2:end), (sum(meanc_t_store(:,2:end).*income_shares)-sum(meanc_t_store_pre(:,2:end).*income_shares))/liqShock, 'LineWidth', 5, 'color', colorB1); hold on;
        load(['output/', benchmark_PB_P, '_HP_shock_25']);
        meanc_t_store_pre = meanc_t_store;
        load(['output/', benchmark_PB_P, '_FP_liquid_HP_shock_25']);
        plot(time_path(2:end), (sum(meanc_t_store(:,2:end).*income_shares)-sum(meanc_t_store_pre(:,2:end).*income_shares))/liqShock, 'LineWidth', 5, 'color', colorPB); hold on;
        
        load(['output/', benchmark_exp]);
        c_preshock = sum(meanc_t_store(:,1).*income_shares);
        load(['output/', benchmark_exp, '_FP_liquid']);
        tmp1=plot(time_path(2:end), (sum(meanc_t_store(:,2:end).*income_shares)-c_preshock)/liqShock, 'LineWidth', 3, 'color', colorB1); hold on;
        tmp1.Color(4) = 0.3;
        load(['output/', benchmark_PB_P]);
        c_preshock = sum(meanc_t_store(:,1).*income_shares);
        load(['output/', benchmark_PB_P, '_FP_liquid']);
        tmp1=plot(time_path(2:end), (sum(meanc_t_store(:,2:end).*income_shares)-c_preshock)/liqShock, 'LineWidth', 3, 'color', colorPB); hold on;
        tmp1.Color(4) = 0.3;
        xlim([0,3]); ylim([0,2]); grid on;
        set(gca,'FontSize', 14);
        xlabel('Years','FontSize', 20, 'Interpreter', 'Latex'); 
        ylabel('Consumption IRF','FontSize', 20, 'Interpreter', 'Latex');
        legend('Exponential','Present Bias','FontSize', 20, 'Interpreter', 'Latex');
        set(gcf,'PaperPosition',[0 0 8 6])
    if saveFigs == 1
        saveas(fig, ['figs/Policy_Functions/AFigure8_1.png'])
    end
    
    fig = figure; hold on;     
        load(['output/', benchmark_exp, '_HP_shock_plus25']);
        meanc_t_store_pre = meanc_t_store;
        load(['output/', benchmark_exp, '_FP_liquid_HP_shock_plus25']);
        plot(time_path(2:end), (sum(meanc_t_store(:,2:end).*income_shares)-sum(meanc_t_store_pre(:,2:end).*income_shares))/liqShock, 'LineWidth', 5, 'color', colorB1); hold on;
        load(['output/', benchmark_PB_P, '_HP_shock_plus25']);
        meanc_t_store_pre = meanc_t_store;
        load(['output/', benchmark_PB_P, '_FP_liquid_HP_shock_plus25']);
        plot(time_path(2:end), (sum(meanc_t_store(:,2:end).*income_shares)-sum(meanc_t_store_pre(:,2:end).*income_shares))/liqShock, 'LineWidth', 5, 'color', colorPB); hold on;
        
        load(['output/', benchmark_exp]);
        c_preshock = sum(meanc_t_store(:,1).*income_shares);
        load(['output/', benchmark_exp, '_FP_liquid']);
        tmp1=plot(time_path(2:end), (sum(meanc_t_store(:,2:end).*income_shares)-c_preshock)/liqShock, 'LineWidth', 3, 'color', colorB1); hold on;
        tmp1.Color(4) = 0.3;
        load(['output/', benchmark_PB_P]);
        c_preshock = sum(meanc_t_store(:,1).*income_shares);
        load(['output/', benchmark_PB_P, '_FP_liquid']);
        tmp1=plot(time_path(2:end), (sum(meanc_t_store(:,2:end).*income_shares)-c_preshock)/liqShock, 'LineWidth', 3, 'color', colorPB); hold on;
        tmp1.Color(4) = 0.3;
        xlim([0,3]); ylim([0,2]); grid on;
        set(gca,'FontSize', 14);
        xlabel('Years','FontSize', 20, 'Interpreter', 'Latex'); 
        ylabel('Consumption IRF','FontSize', 20, 'Interpreter', 'Latex');
        legend('Exponential','Present Bias','FontSize', 20, 'Interpreter', 'Latex');
        set(gcf,'PaperPosition',[0 0 8 6])
    if saveFigs == 1
        saveas(fig, ['figs/Policy_Functions/AFigure8_2.png'])
    end
    
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    
    %Monetary Policy with Income Shocks (Appendix Figure 9)
    load(['output/', benchmark_exp]);
    load(['output/', benchmark_exp, '_Income_shock_5']);
    for tind = 1:length(time_path)
        if tind == 1
            g = g_stationary;
        else
            g = g_t_store(:,:,:,:,tind-1);
        end
        income_shares_t(:,tind) = sum(squeeze(sum(sum(g(:,:,:,:)*da*db))),2);
    end
    fig = figure; hold on;
        load(['output/', benchmark_exp, '_MP_Income_shock_5']);
        meanc_t_store_shock = meanc_t_store;
        load(['output/', benchmark_exp, '_Income_shock_5']);
        plot(time_path(2:end), 100*(sum(meanc_t_store_shock(:,2:end).*income_shares_t(:,2:end))./sum(meanc_t_store(:,2:end).*income_shares_t(:,2:end)) - 1), 'LineWidth', 5, 'color', colorB1); hold on;
        load(['output/', benchmark_PB_P, '_MP_Income_shock_5']);
        meanc_t_store_shock = meanc_t_store;
        load(['output/', benchmark_PB_P, '_Income_shock_5']);
        plot(time_path(2:end), 100*(sum(meanc_t_store_shock(:,2:end).*income_shares_t(:,2:end))./sum(meanc_t_store(:,2:end).*income_shares_t(:,2:end)) - 1), 'LineWidth', 5, 'color', colorPB); hold on;

        load(['output/', benchmark_exp]);
        tmp1=plot(time_path(2:end), 100*(sum(meanc_t_store(:,2:end).*income_shares)./sum(meanc_t_store(:,1).*income_shares) - 1), 'LineWidth', 3, 'color', colorB1); hold on;
        tmp1.Color(4) = 0.3;
        load(['output/', benchmark_PB_P]);
        tmp1=plot(time_path(2:end), 100*(sum(meanc_t_store(:,2:end).*income_shares)./sum(meanc_t_store(:,1).*income_shares) - 1), 'LineWidth', 3, 'color', colorPB); hold on;
        tmp1.Color(4) = 0.3;
        legend('Exponential','Present Bias','FontSize', 20, 'Interpreter', 'Latex');
        xlim([0,3]); ylim([0, 7]); grid on;
        set(gca,'FontSize', 14);
        xlabel('Years','FontSize', 20, 'Interpreter', 'Latex'); ylabel('Consumption Elasticity','FontSize', 20, 'Interpreter', 'Latex');
        set(gcf,'PaperPosition',[0 0 8 6])
    if saveFigs == 1
        saveas(fig, ['figs/Policy_Functions/AFigure9_1.png'])
    end 
    
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    
    %Fiscal Policy with Income Shocks (Appendix Figure 9)
    load(['output/', benchmark_exp]);
    load(['output/', benchmark_exp, '_FP_liquid_Income_shock_5']);
    for tind = 1:length(time_path)
        if tind == 1
            g = g_stationary;
        else
            g = g_t_store(:,:,:,:,tind-1);
        end
        income_shares_t(:,tind) = sum(squeeze(sum(sum(g(:,:,:,:)*da*db))),2);
    end
    fig = figure; hold on;     
        load(['output/', benchmark_exp, '_FP_liquid_Income_shock_5']);
        meanc_t_store_shock = meanc_t_store;
        load(['output/', benchmark_exp, '_Income_shock_5']);
        plot(time_path(2:end), (sum(meanc_t_store_shock(:,2:end).*income_shares_t(:,2:end)) - sum(meanc_t_store(:,2:end).*income_shares_t(:,2:end)))/liqShock, 'LineWidth', 5, 'color', colorB1); hold on;
        load(['output/', benchmark_PB_P, '_FP_liquid_Income_shock_5']);
        meanc_t_store_shock = meanc_t_store;
        load(['output/', benchmark_PB_P, '_Income_shock_5']);
        plot(time_path(2:end), (sum(meanc_t_store_shock(:,2:end).*income_shares_t(:,2:end)) - sum(meanc_t_store(:,2:end).*income_shares_t(:,2:end)))/liqShock, 'LineWidth', 5, 'color', colorPB); hold on;

        load(['output/', benchmark_exp]);
        c_preshock = sum(meanc_t_store(:,1).*income_shares);
        load(['output/', benchmark_exp, '_FP_liquid']);
        tmp1=plot(time_path(2:end), (sum(meanc_t_store(:,2:end).*income_shares)-c_preshock)/liqShock, 'LineWidth', 3, 'color', colorB1); hold on;
        tmp1.Color(4) = 0.3;
        load(['output/', benchmark_PB_P]);
        c_preshock = sum(meanc_t_store(:,1).*income_shares);
        load(['output/', benchmark_PB_P, '_FP_liquid']);
        tmp1=plot(time_path(2:end), (sum(meanc_t_store(:,2:end).*income_shares)-c_preshock)/liqShock, 'LineWidth', 3, 'color', colorPB); hold on;
        tmp1.Color(4) = 0.3;
        xlim([0,3]); ylim([0,2]); grid on;
        set(gca,'FontSize', 14);
        xlabel('Years','FontSize', 20, 'Interpreter', 'Latex'); 
        ylabel('Consumption IRF','FontSize', 20, 'Interpreter', 'Latex');
        legend('Exponential','Present Bias','FontSize', 20, 'Interpreter', 'Latex');
        set(gcf,'PaperPosition',[0 0 8 6])
    if saveFigs == 1
        saveas(fig, ['figs/Policy_Functions/AFigure9_2.png'])
    end
    
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    
    %Monetary Policy with Procrastination Reduction (Appendix Figure 18)
    fig = figure; hold on;
        load(['output/', benchmark_PB_P]);
        plot(time_path(2:end), 100*(sum(meanc_t_store(:,2:end).*income_shares)./sum(meanc_t_store(:,1).*income_shares) - 1), 'LineWidth', 5, 'color', colorPB); hold on;
        load(['output/', benchmark_PB_P]);
        la_slowrefi = -2*log(0.5);
        simulateForward;
        plot(time_path(2:2:end), 100*(sum(meanc_t_store(:,2:2:end).*income_shares)./sum(meanc_t_store(:,1).*income_shares) - 1),  '+', 'color', colorPB);
        load(['output/', benchmark_PB_P]);
        slow_refi = 0;
        simulateForward;
        plot(time_path(2:2:end), 100*(sum(meanc_t_store(:,2:2:end).*income_shares)./sum(meanc_t_store(:,1).*income_shares) - 1), '*', 'color', colorPB);             
    load(['output/', benchmark_PB_P, '_ARM']);
        plot(time_path(2:end), 100*(sum(meanc_t_store(:,2:end).*income_shares)./sum(meanc_t_store(:,1).*income_shares) - 1), 'LineWidth', 3, 'color', colorPB, 'linestyle', ':'); hold on;
    xlim([0,3]); ylim([0, 7]); grid on;
    set(gca,'FontSize', 14);
    xlabel('Years','FontSize', 20, 'Interpreter', 'Latex'); ylabel('Consumption Elasticity','FontSize', 20, 'Interpreter', 'Latex');
    legend('Present Bias (FRM)', 'Present Bias (FRM, $\frac{1}{2}$ Proc.)', 'Present Bias (FRM, No Proc.)', 'Present Bias (ARM)', 'FontSize', 20, 'Interpreter', 'Latex', 'Location', 'Southeast');
    set(gcf,'PaperPosition',[0 0 8 6])
    if saveFigs == 1
        saveas(fig, ['figs/Policy_Functions/AFigure18.png'])
    end

    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    
    
%--------------------------------------------------------------------------
%% Figures -- Intermediate Cases
    
    %Monetary Policy, Beta = 1 (Appendix Figure 10)
    fig = figure; hold on;
        load(['output/', benchmark_exp]);
        plot(time_path(2:end), 100*(sum(meanc_t_store(:,2:end).*income_shares)./sum(meanc_t_store(:,1).*income_shares) - 1), 'LineWidth', 5, 'color', colorB1); hold on;
        
        load(['output/', exp_rational_slow]);
        plot(time_path(2:end), 100*(sum(meanc_t_store(:,2:end).*income_shares)./sum(meanc_t_store(:,1).*income_shares) - 1), 'LineWidth', 5, 'color', [colorB1 0.3], 'linestyle', '--'); hold on;
        
        load(['output/', exp_proc]);
        plot(time_path(2:end), 100*(sum(meanc_t_store(:,2:end).*income_shares)./sum(meanc_t_store(:,1).*income_shares) - 1), 'LineWidth', 5, 'color', [colorB1 0.3], 'linestyle', ':'); hold on;
        
        load(['output/', benchmark_PB_P]);
        plot(time_path(2:end), 100*(sum(meanc_t_store(:,2:end).*income_shares)./sum(meanc_t_store(:,1).*income_shares) - 1), 'LineWidth', 5, 'color', colorPB, 'linestyle', ':'); hold on;
        
        legend('Benchmark Exponential (No Inertia)','Exponential (Rational Inertia)','Exponential (Procrastination)', 'Benchmark Present Bias (Procrastination)', 'FontSize', 20, 'Interpreter', 'Latex');
        xlim([0,3]); ylim([0, 7]); grid on;
        set(gca,'FontSize', 14);
        xlabel('Years','FontSize', 20, 'Interpreter', 'Latex'); ylabel('Consumption Elasticity','FontSize', 20, 'Interpreter', 'Latex');
        set(gcf,'PaperPosition',[0 0 8 6])
    if saveFigs == 1
        saveas(fig, ['figs/Policy_Functions/AFigure10_1.png'])
    end
    
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    
    %Monetary Policy, Beta < 1 (Appendix Figure 11)
    fig = figure; hold on;        
        load(['output/', PB_noproc]);
        plot(time_path(2:end), 100*(sum(meanc_t_store(:,2:end).*income_shares)./sum(meanc_t_store(:,1).*income_shares) - 1), 'LineWidth', 5, 'color', [colorPB 0.3]); hold on;
        
        load(['output/', PB_rational_slow]);
        plot(time_path(2:end), 100*(sum(meanc_t_store(:,2:end).*income_shares)./sum(meanc_t_store(:,1).*income_shares) - 1), 'LineWidth', 5, 'color', [colorPB 0.3], 'linestyle', '--'); hold on;
        
        load(['output/', benchmark_PB_P]);
        plot(time_path(2:end), 100*(sum(meanc_t_store(:,2:end).*income_shares)./sum(meanc_t_store(:,1).*income_shares) - 1), 'LineWidth', 5, 'color', colorPB, 'linestyle', ':'); hold on;

        load(['output/', benchmark_exp]);
        plot(time_path(2:end), 100*(sum(meanc_t_store(:,2:end).*income_shares)./sum(meanc_t_store(:,1).*income_shares) - 1), 'LineWidth', 5, 'color', colorB1); hold on;

        legend('Present Bias (No Inertia)','Present Bias (Rational Inertia)','Benchmark Present Bias (Procrastination)', 'Benchmark Exponential (No Inertia)', 'FontSize', 20, 'Interpreter', 'Latex');
        xlim([0,3]); ylim([0, 7]); grid on;
        set(gca,'FontSize', 14);
        xlabel('Years','FontSize', 20, 'Interpreter', 'Latex'); ylabel('Consumption Elasticity','FontSize', 20, 'Interpreter', 'Latex');
        set(gcf,'PaperPosition',[0 0 8 6])
    if saveFigs == 1
        saveas(fig, ['figs/Policy_Functions/AFigure11_1.png'])
    end
    
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    
    %Fiscal Policy, Beta = 1 (Appendix Figure 10)
    fig = figure; hold on; 
        load(['output/', benchmark_exp]);
        c_preshock = sum(meanc_t_store(:,1).*income_shares);
        load(['output/', benchmark_exp, '_FP_liquid']);
        plot(time_path(2:end), (sum(meanc_t_store(:,2:end).*income_shares)-c_preshock)/liqShock, 'LineWidth', 5, 'color', colorB1); hold on;
        
        load(['output/', exp_rational_slow]);
        c_preshock = sum(meanc_t_store(:,1).*income_shares);
        load(['output/', exp_rational_slow, '_FP_liquid']);
        plot(time_path(2:end), (sum(meanc_t_store(:,2:end).*income_shares)-c_preshock)/liqShock, 'LineWidth', 5, 'color', [colorB1 0.3], 'linestyle', '--'); hold on;
        
        load(['output/', exp_proc]);
        c_preshock = sum(meanc_t_store(:,1).*income_shares);
        load(['output/', exp_proc, '_FP_liquid']);
        plot(time_path(2:end), (sum(meanc_t_store(:,2:end).*income_shares)-c_preshock)/liqShock, 'LineWidth', 5, 'color', [colorB1 0.3], 'linestyle', ':'); hold on;
         
        load(['output/', benchmark_PB_P]);
        c_preshock = sum(meanc_t_store(:,1).*income_shares);
        load(['output/', benchmark_PB_P, '_FP_liquid']);
        plot(time_path(2:end), (sum(meanc_t_store(:,2:end).*income_shares)-c_preshock)/liqShock, 'LineWidth', 5, 'color', colorPB, 'linestyle', ':'); hold on;

        xlim([0,3]); ylim([0,2]); grid on;
        set(gca,'FontSize', 14);
        xlabel('Years','FontSize', 20, 'Interpreter', 'Latex'); 
        ylabel('Consumption IRF','FontSize', 20, 'Interpreter', 'Latex');
        legend('Benchmark Exponential (No Inertia)','Exponential (Rational Inertia)','Exponential (Procrastination)', 'Benchmark Present Bias (Procrastination)', 'FontSize', 20, 'Interpreter', 'Latex');
        set(gcf,'PaperPosition',[0 0 8 6])
    if saveFigs == 1
        saveas(fig, ['figs/Policy_Functions/AFigure10_2.png']);
    end
    
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    
    %Fiscal Policy, Beta < 1 (Appendix Figure 11)
    fig = figure; hold on; 
        load(['output/', PB_noproc]);
        c_preshock = sum(meanc_t_store(:,1).*income_shares);
        load(['output/', PB_noproc, '_FP_liquid']);
        plot(time_path(2:end), (sum(meanc_t_store(:,2:end).*income_shares)-c_preshock)/liqShock, 'LineWidth', 5, 'color', [colorPB 0.3]); hold on;
        
        load(['output/', PB_rational_slow]);
        c_preshock = sum(meanc_t_store(:,1).*income_shares);
        load(['output/', PB_rational_slow, '_FP_liquid']);
        plot(time_path(2:end), (sum(meanc_t_store(:,2:end).*income_shares)-c_preshock)/liqShock, 'LineWidth', 5, 'color', [colorPB 0.3], 'linestyle', '--'); hold on;
        
        load(['output/', benchmark_PB_P]);
        c_preshock = sum(meanc_t_store(:,1).*income_shares);
        load(['output/', benchmark_PB_P, '_FP_liquid']);
        plot(time_path(2:end), (sum(meanc_t_store(:,2:end).*income_shares)-c_preshock)/liqShock, 'LineWidth', 5, 'color', colorPB, 'linestyle', ':'); hold on;
         
        load(['output/', benchmark_exp]);
        c_preshock = sum(meanc_t_store(:,1).*income_shares);
        load(['output/', benchmark_exp, '_FP_liquid']);
        plot(time_path(2:end), (sum(meanc_t_store(:,2:end).*income_shares)-c_preshock)/liqShock, 'LineWidth', 5, 'color', colorB1); hold on;

        xlim([0,3]); ylim([0,2]); grid on;
        set(gca,'FontSize', 14);
        xlabel('Years','FontSize', 20, 'Interpreter', 'Latex'); 
        ylabel('Consumption IRF','FontSize', 20, 'Interpreter', 'Latex');
        legend('Present Bias (No Inertia)','Present Bias (Rational Inertia)','Benchmark Present Bias (Procrastination)', 'Benchmark Exponential (No Inertia)', 'FontSize', 20, 'Interpreter', 'Latex');
        set(gcf,'PaperPosition',[0 0 8 6])
    if saveFigs == 1
        saveas(fig, ['figs/Policy_Functions/AFigure11_2.png']);
    end
    
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    